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1. Introduction. 

The characterization and theoretical understanding of non-equilibrium phenomena 
forms one of the greatest challenges of current statistical physics [HO [3]. In contrast to 
systems in thermal equilibrium, systems far from equilibrium carry nontrivial fluxes of 
physical quantities such as particles or energy. These fluxes are induced and maintained 
by coupling the system to multiple reservoirs, acting as sources and sinks (of particles or 
energy, say) for the system. Since systems of this type abound in biology, chemistry, and 
engineering, any progress in this area is likely to have significant impact, well beyond 
condensed matter physics. 

Over the past years, the remarkable richness of non-equilibrium physics has been 
illustrated beautifully, through the study of simple models. Yet, a satisfactory and 
comprehensive theoretical framework is still missing, even for the arguably simplest 
generalizations of equilibrium, namely, non-equilibrium stationary states (NESS). At 
the source of the difficulties lie several features: First, the evolution of many interesting 
statistical systems is described by a set of transition rates without the basis of a known, 
underlying Hamiltonian. Familiar examples span a wide range, e.g., chemical reactions, 
cell functions, population dynamics, social networks, and financial markets. Second, 
when (and if) such system reaches a stationary state, it will be a NESS in general. 
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For such states, there is no a priori knowledge of, or a good working hypothesis for, 
the configurational probability distribution. By contrast, for systems known to reach 
thermal equilibrium, the fundamental hypothesis is highly successful in providing the 
stationary distribution, regardless of the details of the dynamical evolution. Finally, 
seemingly minor modifications of the evolutionary rules (i.e., the system dynamics) 
often lead to dramatic changes in macroscopic properties, suggesting that very different 
NESS's can be associated with just slightly different dynamics. 

To appreciate how profoundly NESS differ from equilibrium systems in this regard, 
let us briefly recall the Gibbs-Boltzmann framework for systems in thermal equilibrium 
(with a single reservoir). For simplicity, we use the language of the canonical ensemble 
here The required inputs are: a set of configurations C±, C 2 , ... of the system 
(also known as microstates) and an expression for the internal energy TC(C) associated 
with each configuration C. Then, the statistical weight, P eq (C), for a system in thermal 
equilibrium with a reservoir is well established, in terms of T~C{C) and a simple parameter 
- temperature T - associated with the reservoir. All macroscopic observables now follow 
as configurational averages, (A) = J2c A{C)P eq {C). In particular, a system coupled to 
a heat bath at inverse temperature (3 = l/ksT is governed by the familiar Boltzmann 
distribution, P eq (C) = Z~ x exp[-pH(C)]. 

Of course, all real systems, whether in equilibrium or not, are fundamentally 
dynamic in nature, continuously undergoing transitions from one configuration to 
another. In that context, stationary distributions, including P eq (C), emerge as the 
long-time limit of a time-dependent distribution, P(C; t). For its time evolution, we will 
assume that it obeys a master equation, with a given set of transition rates between 
the configurations. From a modelling perspective, it is therefore natural to ask: What 
choices of dynamic transition rates will lead to a desired stationary distribution? and: 
What modifications of these rates will leave this distribution invariant? For equilibrium 
systems, the answer is well known and can be traced to the property of detailed balance, 
which is related to microscopic reversibility. Specifically, if the rates satisfy detailed 
balance, the net probability current between any pair of configurations vanishes in the 
steady state. As a result, P eq {C) = lim^oo P(C; t) can be expressed in terms of ratios 
of these rates, and the long-time limit remains invariant under any modification of 
the dynamics which preserves these ratios. Indeed, Monte Carlo simulation studies of 
equilibrium systems rely heavily on this property. Since the full configurational sum 

is computationally inaccessible, an importance sampling of configuration space is 
achieved dynamically, by constraining the transition rates such that P(C; t) eventually 
approaches the desired equilibrium distribution. 

In this article, we begin by addressing two fundamental questions. First, is there 
a general procedure to find the stationary solution(s) of a master equation even if 
the transition rates violate detailed balance? Second, can we specify the class of 
transformations of the rates which leave this stationary distribution invariant? The 
answer to the first question is not new, expressing the solution in terms of directed trees. 
Since it does not seem to be widely known, we will provide a brief review here, restricting 
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ourselves to systems with (i) a finite (but arbitrary) number of configurations, (ii) time- 
independent transition rates, and (Hi) ergodicity. In this framework, we discuss the 
importance of irreversible loops (associated with the transition rates) in configuration 
space, the role they play in the non-trivial probability currents, and how equilibrium 
distributions are recovered in case all loops are reversible. To answer the second question, 
we proposed j5] a general classification of NESS in terms of the stationary configurational 
probabilities P*(C) and the (stationary) probability currents connecting them, K*(C',C) 
- denoted in short by {P*,K*}. A given set {P*,K*} defines a particular NESS, and 
allows us to compute all physical currents - mass, energy, etc - for this state. We then 
discuss the set of transformations of the transition rates which leave {P*, K*} invariant. 
In other words, we describe how to generalize the "detailed balance condition" to NESS. 

In the reminder of this article, we explore some consequences of our postulate and 
present several specific examples, to illustrate the very general and formal framework 
proposed. In a concluding section, we provide a summary and outlook, including some 
remarks on generalizations to discrete time and continuous configuration space. The 
Appendices are devoted to technical details. 

2. The master equation and its associated steady state 

We first establish the mathematical framework for our analysis. Consider a general 
continuous-time dynamics, with a discrete, finite configuration space. We assume that 
every configuration can be reached from every other one, so that the system is ergodic. 
Labelling the configurations in some arbitrary manner as C\, C2, •-, Ctv, we are interested 
in P (Cf, t), the probability to find the system in configuration at time t. Its evolution 
is governed by a set of transition rates w (Cj —>■ C»), for the system to change from 
configuration Cj to configuration Cj, per unit time. To simplify notation, let Pi(t) stand 
for P (Ci,t) and w (Cj — > Cj) be denoted by [Bj. All w\ are real, non-negative, and 
assumed to be time- independent. In general, "forward" and "backward" rates differ, 
i.e., w\ 7^ w l j. If one of them vanishes while the other remains nonzero, we will call 
the corresponding transition uni-directional. In terms of these to's, the master equation 
simply expresses the rate of change of Pi(t) as the system makes transitions in and out 



of Ci. 




(1) 



This is often more conveniently written in matrix form 



N 




(2) 



where we have introduced the N x N matrix W via 




(3) 
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We note that W- is a stochastic matrix (in the continuous-time formalism) since it 
satisfies (i) W\ > for all i ^ j, and (ii) Y2i W{ = for all j. The second condition 
ensures that, once normalized, Pi(t) remains so for all subsequent times. 

The right hand side of Equation ([1]) can be expressed in terms of the (net) probability 
currents Kf, from configuration Cj into configuration Ci, 

K{ (t) = wiP.it) - w'Piit) (4) 

so that the master equation simply states the conservation of probability: d t Pi(t) = 

Since the dynamics is ergodic, Equation jT]) has a unique stationary solution, 
P* = lim^oo Pi(t), independent of the initial conditions. This implies that P* is a right 
eigenvector of W with eigenvalue zero. This eigenvalue is non-degenerate so that P* 
spans the null space of the matrix W- . The associated stationary currents are denoted 
by K*\. They satisfy the equality = 0> f° r au h i- e -> the total probability 

current into any given configuration vanishes. 

In general, P* depends on the chosen rates, w\. When simulating systems in thermal 
equilibrium, the challenge is to specify a set {w?} such that the resulting stationary 
state equals the desired equilibrium distribution, P[ q . A well-established procedure is 
to choose rates which satisfy "detailed balance" (with respect to -Ff 9 ), namely, 

^p; q - ^p- q = o (5) 

for every pair Ci, Cj of configurations. This relation can be viewed as a constraint on the 
set of allowable io's. We see immediately that Equation (jHJ) is equivalent to demanding 
that all individual currents vanish, i.e., {K eq ){ = for all i ^ j. Conversely, if the 
dynamics satisfies detailed balance, even if the steady-state distribution is not explicitly 
known, we can construct P^ q easily from Equation ((Sj). Furthermore, it is now very easy 
to determine whether two different sets of rates, {w{} and {w{}, will generate the same 
P e - q : They do if the equalities wf/Wj = wl/m* hold. 

Equation (jSJ) seems to imply that P^ q must be explicitly known in order to test the 
validity of detailed balance. However, detailed balance is an intrinsic property of the 
dynamics, expressing a deeper statement on microscopic reversibility and requiring no 
information about any specific steady-state distribution P^ q . Known as the Kolmogorov 
criterion [3 El 121 IS] , ^ relies on considering closed loops in configuration space, e.g., C 
= Ci — > Cj — > Ck — > ... — > C n — > Ci. For each loop, we define the product of the 
associated rates in the "forward" direction, II [C] = w l jW 3 k ...w™ , as well as for the 
"reverse" direction: II [C rev ] = wlvjj ... w l n . In terms of these products, the dynamics is 
said to satisfy detailed balance if 

n [£] = n [c rev ] (6) 

for all loops. This condition implies the path-independence of the ratio of the associated 
products of the rates along any path which goes from Ci to Cj. To be more explicit, 
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consider a path: V = Cj — >■ ... — > C n — > Cj and its "reverse" P rei) = Cj — > C n 
Cfc — > Cj. Then, the ratio 



wl...w? 



W J n ...wf 

assumes the same value as II* [P']/ IT^ [T^eJ where V is any other path from d to Cj. 
Such path independence allows us to define a unique "potential" associated with each 
Cf. 

$ io = in(n°/n l ) , (8) 

where C Q is some arbitrary reference configuration. As we will demonstrate below, it is 
straightforward to show that the stationary distribution is then given by 

P* oc exp [$ io ] . (9) 

The relationship between this P* and the familiar P[ q is now clear: The potential $ 
is, e.g., just — P {H [Ci] — H[C ]} for the canonical distribution. Further, for systems 
obeying Equation ([6]), this approach allows us to define an effective Hamiltonian "H e ff 
along with an effective (inverse) temperature /3 e //, m case such concepts are beneficial 
for the problem at hand. 

A key signature of non-equilibrium steady states is that the underlying rates do 
not satisfy microscopic reversibility and Equation ([6]) is violated. As a result, the ratios 
Equation ([7j) will be depend explicitly on the path V . This leads Lebowitz and Spohn 
to define [10] an action associated with each V as the log of the ratio ((71). In a more 
specific arena, the dynamic functional [11] is its counterpart in typical statistical field 
theories. 

For a NESS, the right hand side of Equation (jSJ) is generically nonzero, given by the 
nontrivial stationary probability currents: As a result, the construction of the stationary 
probability distribution - while still possible - requires a bit more effort. Though 
established some time ago [121 131 E] > the graph-theoretical approach, similar to those 
originally designed for electric networks [15], appears not to be widely known [9]. For 
the reader's convenience, we briefly review the method here. First, we represent each 
configuration by a labelled vertex (e.g., i for d). Next, we draw all distinct spanning 
trees (i.e. all distinct graphs containing all vertices and exactly one single undirected 
edge between each pair so that no loops are formed) t a , a = 1, 2, .., M, with iV vertices. 
The total number of such trees, M, is given by N N ~ 2 , according to Cayley's theorem 
|16j . Next, we select a specific vertex, e.g., i, and draw an arrow, directed towards i, on 
every edge. The resulting directed tree will be denoted by t a ^. Note that, for given i, 
there is exactly one directed tree for each undirected tree. Next, we associate a factor 
Wj with an edge directed from vertex k to j. Finally, the numerical value U(t a ^) of the 
directed tree t a u\ is defined to be the product of the (N — 1) rates w appearing in that 
particular tree. If one of the rates vanishes, we simply assign the value U{t a ^) = to 
the associated tree. We illustrate the procedure in Figure [H for finding P x * in a iV = 4 
case. While we show all 16 directed trees, we just give two of the C/'s as examples: 
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Figure 1. All 16 directed trees contributing to the graphical representation of P*, for 
a simple model with N = 4. 

U{ti(i)) = wfwfwl (top right diagram in Figured]) and U(ti 3 ^) = w\w\w\ (bottom left 
diagram in [I]). The stationary solution of Equation ([1]) is then given by 

P* = Z-^Ufav) (10) 

a 

with the normalization factor Z defined as 

N 

z = J2Y. u m ( n ) 

i=l a 

Now that we have a representation for P*, let us consider the probability currents, 
defined in Equation (TJJ above. Specializing to the stationary case, we write the net 
current from Cj into Cj in the form: 

K*{ = w(P* - w)P* = Z- 1 HU(t aU) ) - wp(t a(i) )] (12) 

a 

Focusing on the expression within [...], we note that, for a specific a, the trees t a ^ 
and differ only in the directed edges that connect vertices i and j. In Figure [21 
we illustrate this statement with a tree that has k\, hg as the vertices between i to 
j. In Figure [2k,, t a ^ involves a directed tree used for U{t a ^), with directed edges 
being blue. Similarly, in Figure [2b, we have the same tree being used for U(t a ^), the 
only difference being the edges between i and j which are now reversed. Considering 
just this segment of the two trees, we may write the associated products of the w's as 
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(a) (b) 

Figure 2. An illustration of the forward and backward loops, appearing in the sum 
on the right hand side of Equation p^|) . See text for details. 



n*(t Q (j)) = w t kl w , j£ 2 ...Wj t and U.{(t a ^) = w J k ...w^ 1 , respectively. Meanwhile, the rest of 
both trees (the side branches "coming into the circle" in Figure [2]) are identical. Thus, 
denoting the products of the rates associated with these side branches as R, we may 
write 

Combining the side branches with the path between % and j, we arrive at 

U(t a{j) ) =TV J (t aU) )R(t a{j) ) and U{t a{i) ) =TV l {t a{l) )R{t a{i) ) (13) 

Returning to the steady-state current, Equation (fT2l) . we see that the additional factor 
of w in each term can be regarded as adding an extra edge (red), so that each tree is 
converted into a graph with a single loop (the "circles" in Figure [2]) . It is natural to 
label these directed loops as C a (j) and C a (i) respectively. Of course, both refer to the 
same loop, just traversed in opposite directions. Quantitatively, we have, for each t a 

™iU(t a(j) ) - wp(t a(i) ) = {wiwi 1 wl]...wf - w i j w 3 kt ...w$ 1 } R(t a(j) ) . 

Now, the terms in {...} can be associated with C a (j\ and C a n), so that a compact 
expression for the current is 

K*\ = Z- 1 [n(£oy)) - n(£ a(i) )] R(t a[i) ) . (14) 

a 

From here, it is immediately obvious that detailed balance, manifested in reversible 
loops, Equation (jSJ), will lead to all K*\ being zero, and nontrivial steady-state currents 
can only emerge from rates which violate detailed balance. For completeness, we should 
remark that it is conceivable for accidental cancellations in the sum over trees to occur 
so that some K*'s may vanish even though irreversible loops are present. However, K* 
should be nonzero, generically. 

Let us point out briefly that there is an alternate representation of the stationary 
distribution [9j, in terms of the co-factors of W- (which we denote by Cj here). Of course, 
a well known expression of the determinant leads us to det W = . W-Cj for any i. 
Moreover, for a stochastic matrix W- , Cj is independent of i. Thus, we may write Cj 
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in place of Cj. Finally, since W has a zero eigenvalue, we have = det W = £\ W-Cj. 
Thanks to the uniqueness of the stationary solution, this equality implies Cj oc P*. 

Before proceeding to the next section, it may be instructive to see how the complex 
expression for P*, Equation ffTUl) . reduces directly to the familiar form for systems in 
equilibrium. To start, we choose an arbitrary reference configuration C Q . Of course, we 
will rely on the path- independent properties of U°/W . Consider the ratio 

a 

and exploit the factorization, Equation (fTBl) . Then, the numerator of the above reads 
E U M = E n°(t a(j) )R(t a(j) ) , (15) 

a a 

where iK(t a (j)) now stands for the product of the rates associated with the path from C Q 
to Cj in t a u-\ and R(t a u-\), with the rest of the tree ("side branches"). Of course, we have 
a similar decomposition for the denominator. Next, we invoke the path independence 
of the il's (Equation (jSj)) and write 

Yl°(t a(j) ) = Eg(fa(o))e** • 

Let us emphasize that $j G does noi depend on the path so that it is also independent 
of t a . Consequently, Equation ffl5|) becomes 

a a 

But, as noted above, the side branches of t a ^ and t a (o) are identical, so that we can 
replace R{t a ^) by R(t a ^). The sum on the right is easily recognized as J2 a U(ta(o)) 
and we arrive at 

E^0))=e*-E f/ (^)) ' 

a a 

and the desired result for systems in thermal equilibrium: 

i^ = P;exp[$ lo ] . (16) 



E^w) 

a 



3. A postulate for the complete characterization of NESS 

In much of the literature on NESS, the microscopic stationary distribution P* tends to 
be the central focus, along with those macroscopic properties which can - at least in 
principle - be derived from it, such as order parameters and correlation functions. In 
this sense, the investigations follow standard routes for systems in thermal equilibrium. 
Indeed, many studies emphasize the similarity between the P* of a NESS and the 
P eq of an equilibrium system. At the same time, many complementary studies focus 
on quantities absent from equilibrium systems, e.g., particle currents or energy fluxes 
through the system. Unlike P*, however, these are macroscopic averages, and thus 
located at the same scale as order parameters and correlation functions. It is natural 
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to ask whether there is a non-equilibrium quantity, or concept, which would provide a 
microscopic basis for these fluxes. Here, we propose that the microscopic distribution 
of probability currents, K*, be brought to stage center, as an indispensable partner for 
the microscopic probability distribution P*. Clearly, the particle/energy fluxes can be 
computed from K* as averages, on a par with (•) = £ • P*. Needless to say, in the 
arena of equilibrium statistical physics, all K*'s vanish and therefore play the role of an 
invisible partner. 

To showcase the essential nature of K*, let us consider a completely trivial example, 
namely, a single particle hopping randomly between neighboring sites on a ring. As long 
as the rates are spatially uniform, the stationary probability distribution is trivially flat: 
P* oc 1, regardless of whether the hopping rates are symmetric (i.e., identical for left and 
right hops) or biased. Yet, detailed balance is satisfied if and only if they are symmetric, 
and there is an important physical difference between the two cases, namely, the presence 
or absence of a particle current in the stationary state. In more complex cases, other 
types of currents (e.g., energy) may be the key feature. In general, all these fluxes can 
be traced to the violation of detailed balance and the associated microscopic probability 
currents. In other words, without K*, P* alone cannot be a complete description of a 
NESS. 

That K* has been largely ignored in the literature may be traced to the following. In 
all studies of non-equilibrium systems, we begin with a given set of rates {w}, motivated 
by the underlying physics, chemistry, biology, psychology, sociology, etc. The main 
difficulty is to find P*. But, once P* is known, macroscopic average fluxes can be 
computed easily from w and P*, so that it is unnecessary to construct K* explicitly. 
However, if we wish to ask the inverse question, namely: "What rates are needed to 
achieve a particular NESS?" , then K* plays an indispensable role. 

Though seemingly trivial, let us state for completeness that K* alone is also 
inadequate. Specifically, it is possible to have NESS's with the same K* but different 
P*'s. Electromagnetism provides a good analog. In electrostatics, the central focus is 
the charge distribution, being the analog of P*. By definition of electrostatics, there 
are no currents anywhere. Similarly, we are concerned mainly with the currents in 
magnetostatics and tend to ignore the charge distribution. However, it is clear that, in 
general, the charge distribution is not trivial, especially for non-neutral systems. 

Staying with electromagnetism for a moment, Kirchhoff's solution for general 
electric circuits [15J is often quoted as the first graphical solution to the master equation. 
Yet, there are non-trivial differences. Specifically, the currents (iy between nodes i and 
j) are clearly the quantities of interest for Kirchhoff, while there is no trace of the 
charge distribution (p, at node i) in his solution. In contrast, P* is the key quantity 
when solving the master equation. Further, there is no one-to-one mapping from the 
set of electromotive forces and resistances to the set of rates. Further details regarding 
the "duality" of our NESS and the Kirchhoff problem may be found in Appendix A. 

Motivated by these considerations, we propose that {P*,K*}, the distributions 
for probability and probability currents, form a complete and unique description for 
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any stationary state. This classification includes equilibrium systems, characterized by 
{P eq ,0}, as well as NESS's where K* is, by definition, non-zero. 

The inclusion of K* in this characterization motivates the definition of the 
"distance" of a given NESS from equilibrium, namely, 

\K*\ 2 = J2( K 1) 2 ■ ( 17 ) 

hj 

This may provide a quantitative basis for phrases like "near equilibrium" and "far from 
equilibrium." How useful this concept is remains to be explored. In Section fl~2l below, 
we will see (K*^) 2 appearing in a discussion of entropy production. 

More importantly, K* contains all information necessary to determine the average 
fluxes (or currents) associated with physical observables, such as energy or particle 
number density. A microscopic distribution of probability currents, it serves as the 
statistical weight in the computation of physical currents. Explicitly, we write 

where J (Cj, Cj) is associated with a physical observable. Naturally, we expect physically 
meaningful jT's to be antisymmetric under Ci <s> Cj (e.g., a particle current). Taking 
advantage of this fact, we can define (J) in a manner that highlights K* as a 
probabilistic weight, namely, 

(J) = Y,J (19) 

{id} 

where {i,j} denotes a sum taken over only the positive K*s. For the readers' 
convenience, we illustrate these considerations with two examples. 

Many NESS models involve particles hopping from a site (labeled by s) to a 
neighboring site on a lattice, with an excluded volume constraint. A classic example 
is the totally asymmetric exclusion process (TASEP) [T7J [18]. For such models, a 
configuration is uniquely given by the set of occupation numbers: {n s }. Now, one 
quantity of interest might be the average particle current across a particular pair of 
sites (bond). Specifically, for the current from site a to b, we need 

Jbond ({n s } , {n' s }) = n a (1 - n b ) (1 - n' a ) n' b -n' a {\- n' b ) (1 - n a ) n b (20) 

Clearly, for the net current across, say, some surface in a 3-dimensional system, 
we only need to sum the J's associated with bonds which pierce the surface: 
Jsurface {{ n } j { n '}) — ^ ' Jbond- Further, if the entire distribution for this current 
(denoted by p (J)) is of interest, then it is given by 

P^ = l E HJ- Jsurface ({n},{n'}))K*({n},{n'}) . 
M,{«'} 

Let us reiterate the need for K* (as opposed to P*) in finding these physical currents, 
even though P* alone may appear to be the only necessary ingredient for computing 
averages. This paradox can be traced to Equation (j4]), which allows us access to K* from 
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Figure 3. A simple kinetic Ising chain, coupled to two temperature baths. Spins 
coupled to the higher (lower) temperature are indicated by red (blue) dots. In the 
steady state, there are constant energy fluxes, as denoted by the black arrows. 

P*, provided the rates {w} are known. Typical studies of NESS begin with a given set 
of rates and unknown {P*, K*}. Thus, the currents are seemingly unimportant, once P* 
is found. However, as will be discussed in the next section, there are many w's that can 
lead to the same {P*, K*}\ In this sense, the rates are not necessary (though sufficient) 
for finding average fluxes, and some details of their precise form are "irrelevant." 

Another example of a current in a simple NESS is the total energy flux through a 
kinetic Ising model coupled to two thermal baths at different temperatures, as studied 
in, e.g., [19l [20], [2H HI [22J, [23] . For simplicity, we consider a one-dimensional version here, 
as shown in Figure [H Let us denote the set of sites coupled to the two baths by {s} 
and {§}, respectively, and a particular spin configuration by {o~ s }- We are interested in 
the energy flowing from the bath at the higher temperature, T s , to the bath at lower 
temperature, Tg. Given the details of the coupling, energy flows into (out of) the system 
through the {s} ({§}) spins. Using a random sequential dynamics, the energy current 
operator, in units of the nearest-neighbor coupling strength, can be written as 

s 
s 

where S s = Ylk^s (1 + °&°fc) /2 an d a similar 5$ insure that only the spin at site s or 
s is flipped. If the two temperatures are the same and the system is in equilibrium, 
then (J £ ) is trivially zero due to K* = 0. If we use the rates and P* instead, a tedious 
calculation will show that the average of each term in J £ vanishes. 

To summarize, we postulate that {P*,K*} is a complete characterization of a 
general NESS. Compared to an equilibrium state, which is specified by N — 1 quantities 
(P*), a NESS requires O (N 2 ) quantities. Note that the N (N - 1) /2 quantities in K* 
are not independent: There are iV constraints, namely, T^K* 1 - = 0. The significance 
of our characterization is that {P*,K*} gives access to all macroscopic quantities of 
interest, including any currents flowing through, or current loops within, the system, 
without requiring any additional knowledge about the rates. In the following, we explore 
the consequences of our proposal. 
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4. Consequences of the postulate 

Once a NESS is characterized by {P*, K*}, we can explore the possibility of a generalized 
"detailed balance condition," which the transition rates must satisfy in order to 
ensure settling into a given NESS. This freedom of choice has implications for entropy 
production, as well as for the design of optimized computer simulations. 



4-1. Dynamic equivalence classes. 

In many investigations of NESS, one of the most striking features is their sensitivity 
to the details of the transition rates. Major changes of macroscopic properties often 
emerge from seemingly minor modifications of the rates. It is therefore natural to ask 
how one might determine whether two sets of rates will lead to the same, or a different, 
NESS. This issue can be addressed within the framework proposed here. In particular, 
since all properties of a NESS are supposedly given by {P*,K*}, all time-independent 
(static) quantities of physical interest can be computed, without further recourse to the 
dynamics [21] . The analog in equilibrium systems is that all (static) quantities can be 
obtained from P* without detailed knowledge of the rates. 

An alternate phrasing of the above question is: Given a set of u>'s and its 
associated NESS, what are the transformations on the rates which leave this {P*,K*} 
invariant? For the equilibrium case, {P eq , 0}, the answer is provided by the detailed 
balance condition, Equation (J5j): Any set of w's will lead to the desired P eq , provided 
wj/w l j = P^ q /Pj q . Alternatively, this relation can be interpreted as a constraint to 
be placed on the rates if the goal is to reach a specific P eq . In our framework, this 
constraint can now be easily generalized: To arrive at a given {P*,K*} final state, the 
w's must satisfy 

w{p; -w)P* = K*l . (21) 

for all pairs % ^ j. As a simple extension of the equilibrium case (Equation [5]), 
this equation constrains, say, the "backward" rate once an (arbitrary or suitable) 
"forward" rate w\ is selected. In this sense, there is just as much freedom in choosing 
rates to arrive at a given NESS as for systems to reach thermal equilibrium. In the 
following, we will explore alternative expressions for this simple constraint, to provide 
further insight. 

Motivated by the antisymmetric nature of the currents K*j, we define 

W> = W}P* (22) 

and decompose it into its symmetric and antisymmetric parts: 

Wi = S(+4, (23) 

where Sf = (Wf + W]) /2, and A\ = {W( - W]) /2. Then, Equation (f2T|) is a very simple 
constraint, namely, 

4 = \k*{ ■ (24) 
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As a result, for a given {P*,K*}, A is completely determined. In contrast, S is 
essentially free, except for two restrictions. First, the physical rates must be non- 
negative (w > 0), leading to Sf > \Al\, for all i ^ j. Further, probability conservation 
imposes Ej W- = 0, for all j, so that Ej W- = also. Now, Ej A] also vanishes, since 

*Ei(Wf-Wj) = - E* which is zero b y virtue of = °- Therefore, Ej = 

reduces to Ei 5"/ = and we arrive at the restrictions for S: 

S}>l\K*{\ V^j, .S-; - »/ . (25) 
Within these constraints, we may choose arbitrary S's, construct transition rates via 



si + \k*\ 



(p*r , (26) 



and be certain that the final NESS remains unaffected. In this respect, we may regard 
all such S"s - specified by iV (N — 1) /2 parameters - as generating an "equivalence class" 
of dynamic rates associated with the same NESS. 

There is another, perhaps simpler, way to express this freedom of choice. Suppose 
we found P* from a given set of rates w\ (and so, K* is also determined). To construct 
another set of equivalent rates, it is sufficient to add to the w's any set of changes A 
that satisfy 

Aip; = A;/y . (27) 

Note that A may be negative, provided the new rates w\ + A^ are non-negative. 
Of course, this amounts to the same statement as Equation (|2"5"|) . Equation (j2"7|) is 
reminiscent of the "ordinary detailed balance" condition. The distinguishing feature 
here is that the differences between two sets of rates, as opposed to the rates themselves, 
must satisfy "detailed balance with respect to P*." 

Another interesting corollary allows us to determine (some aspects of) P* from two 
different sets of rates w, w' provided we know that both lead to the same NESS. The 
key is to compute the differences: A = w — w'. For all nonvanishing pairs {A], A* }, 
the ratio can be used to construct P*. In this sense, the A's can be thought of as 
'rates that lead to an equilibrium state', satisfying microscopic reversibility, Equation 
06]). Two features distinguish this case from a true equilibrium system: (a) Some A's 
can be negative and (b) there is a unique P* even if the A's are not ergodic (i.e., not 
all configurations connected by A's). 

Lastly, let us present an alternate approach. For equilibrium steady states, instead 
of considering W, it is more convenient to define another matrix, namely, 

wi = wi(p*ip:) 1 ' 2 (28) 

Due to detailed balance, W is symmetric, and the master equation is transformed into a 
standard eigenvalue problem. Apart from a single zero eigenvalue (associated with P eq ), 
all other eigenvalues are negative. Clearly, we can also consider the properties of W in 
the non-equilibrium case. Similar to W, W can be decomposed into a symmetric (S) 
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and a nonzero asymmetric part, A = K*/2. Here, K* is related to the original currents 
K* via 

K*i = K*i(p*p*y 1/2 . (29) 

For the NESS, we can still conclude, thanks to the Perron-Frobenius theorem, that 
there is only a single zero eigenvalue, and that all other eigenvalues have strictly 
negative real parts. However, the eigenvalues need no longer be real and may appear in 
complex conjugate pairs. These correspond to oscillatory relaxation, as in underdamped 
oscillators. 

Similar to our approach above, for a given NESS, A is completely determined while 
S is essentially arbitrary (apart from the minor restrictions imposed by positivity and 
normalizability). The two approaches have complementary advantages: W provides 
immediate information on the probability currents. W has the same spectrum as the 
original W, and therefore reflects the dynamics more clearly. 

We conclude this section by emphasizing that W matrices belonging to an 
equivalence class have, in general, a different set of nonzero eigenvalues. Needless to say, 
the physical interpretation is that different rates generally lead to different relaxation 
rates into the same NESS. Indeed, our hope is that future Monte Carlo simulations 
of the NESS will exploit this freedom to devise more efficient algorithms, in the same 
spirit as, say, the cluster algorithms [25], designed specifically to circumvent critical 
slowing-down in equilibrium systems near critical points. 



4-2. Entropy production 

One of the key signatures of non-equilibrium steady states, recognized many decades 
ago [26j [13l E7] , is entropy production. Using the expression for average fluxes above 
(Equation [T8|) . one possible definition [13] for the total entropy production is: 

S«( t )4j>l( t )> n™. (30) 

i,j 

If the master equation is interpreted in the language of chemical reactions, 
In (W- Pj/WjPfj would be an affinity, or generalized "thermodynamic force" |13j . 
Inserting Equation (TjO) for K\{b), it is immediately apparent that S to t(t) > 0. Further, 
S tot can be written as the sum of two terms, i.e., the entropy production of the "system" 
and of the "medium" , 

i p (+) i ]xri 

Ks(t) ^2^ m)ln pM ' * med{t) E 2^^ (t)ln ^' (31) 

i,j i,3 3 

The former is readily recognized as the time derivative of S sys = — £\ Pi(t) lnPj(t), 
which motivates the term "entropy production of the system" . The latter is attributed 
to the coupling of the system to the external environment in a manner that prevents 
it from reaching equilibrium [13]. Unlike S tot , neither S sys nor S me d are necessarily 
positive. 
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Since we are primarily interested in steady states, we take the infinite time limit. 
For equilibrium systems, all K*{ are identically zero so that both S sys and S me d vanish. 
For a NESS, however, only S* sys vanishes. In general, K* ^ leads to S med = S^ ot > 0. 
The interpretation of these results is clear: In the steady state, the entropy associated 
with the system no longer changes. However, being coupled in an irreversible way to 
the environment, a NESS continues to increase the entropy of its surrounding medium. 
In other words, S med carries information about the precise nature of these couplings, 
encoded in the transition rates. So, even if we insist on having the same NESS (i.e., a 
given {P*,K*}), S med will not be unique. In the following, we explore the connection 
between S med and the transition rates in more detail. 

In the steady state, S£ ot = S med . Recalling the decomposition of WP* into a 
symmetric and an antisymmetric component, Equation ( l23l) . it is more convenient to 
work with S^ ot whence 

i,j J 1 i,j i i 

While the antisymmetric component is fixed by the currents, A = K*/2, the symmetric 
part S can be chosen at will, as long as it satisfies Equation fl25|) . In particular, while 
S^ ot > has to remain valid for any choice of S, it is possible to make the entropy 
production arbitrarily small or infinitely large. In other words, exploiting the freedom 
in the choice of the transition rates, we can select an arbitrary value of S^ ot , without 
affecting the underlying NESS {P*,K*}. Let us provide a few details. 

To achieve an arbitrarily small S£ ot , we only need S 3> A, for every non- vanishing 
element. Expanding Equation ( l32l) to lowest order in A\jS\ = K* 3 i /(2S^), we arrive at 

.2 r 



(kit 



K* 



(33) 



S 

Of course, even though S^ oi can be made arbitrarily small, it still remains strictly 



2 

positive, retaining the NESS signature thanks to (K*fy being strictly positive. Since 
S me d = S toi = characterizes an equilibrium system, minimizing S^ ot (for a given NESS) 
corresponds to selecting rates which are as "equilibrium-like" as possible. 

At the opposite extreme, we can construct rates with "infinite" S med by reducing 
at least one off-diagonal Sf to \ A\\, so that either Sf + A\ or S{ — A\ vanishes (but never 
both). Translating this back into a new matrix of transition rates, via Equation fl26|) . 
the corresponding transition is now uni-directional, in that one of the two directed edges 
between the associated pair of configurations is missing. Of course, it is possible to make 
all rates uni-directional, which may be naturally referred to as "maximally asymmetric." 
Such systems appear frequently in the literature, a particularly familiar example being 
the totally asymmetric exclusion processes (TASEP) [TTIIH]. One clear advantage of 
having maximally asymmetric rates for all edges is that the number of nontrivial trees 
contributing to P* is kept at the absolute minimum. Needless to say, the expression for 
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1 2 




3 4 

Figure 4. All configurations of a TASEP on a one-dimensional lattice with just two 
sites. Particles are denoted by black circles. 

K* also simplifies, to just one term in Equation ffl2]) . and all irreversible loops are also 
uni-directional. 

5. Examples 

The formalism we presented is quite general. It is instructive to provide a series of 
simple examples, to illustrate how it applies in various circumstances. We will start 
with a minute system, a special case of the N = 4 example of El It is closely associated 
with the zero range process (ZRP) and its generalizations which we consider in the 
second subsection. The third subsection is devoted to a kinetic Ising model coupled to 
two thermal baths at different temperatures. Finally, we close with the general class of 
NESS with Gaussian distributions. 

5.1. TASEP with two sites 

The totally asymmetric exclusion process (TASEP ) on a one- dimensional chain with 
open boundaries [T7J Q2] enjoys considerable attention as a non-trivial system with a 
known NESS distribution. The dynamical rules are as follows: Particles may enter 
(leave) a lattice with L sites at the left (right) end with rate a (j3), provided the first 
(last) site is empty (occupied). Within the system, particles may hop to the right, with 
rate 7, provided the target site is empty. Though most of the interesting properties 
appear in the thermodynamic limit, our goal here is to illustrate the considerations of 
the previous sections. For that purpose, it is more helpful to study a small system. To 
ensure that nontrivial loops can exist in configuration space, the smallest 'interesting' 
TASEP is the one with two sites (L = 2) and 4 configurations. Let us label these as 
follows: 1 (4) for both sites being empty (full) and 2 (3) for only the right (left) site 
being occupied (Figure HI). In Figure E^a), we show all the allowed transitions (arrows). 
The rates are a, (3, and 7 for the vertical (blue), horizontal (green), and diagonal (black) 
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Figure 5. For a TASEP with two sites, all non-vanishing wij are indicated by arrows 
in (a), (b) The only directed tree contributing to P*. (c,d) The two directed trees 
contributing to P 3 *. 



arrows. The associated matrix W rates can be written easily: 



W 



( 


—a 


P 





\ 







- (a + (3) 


7 







a 





-7 


p 


\ 





a 








(34) 



Note that it is "maximally asymmetric" , in concordance with every transition being uni- 
directional. As a result, most of the 16 possible trees (cf. Figured]) do not contribute 
to the graph-theoretical representation of the stationary probabilities P*. For example, 
only a single tree remains in the representation of P± (shown in Figure [5(b)), and only 
two trees (Figures [5](c) and [5(d)) remain for P 3 * (see Appendix B for details). The 



nonvanishing trees are easily found: -P 1 * 2 3 4 oc P r yf3,a^,aa(3 + j3(3a, } a r ya.. For clarity, 
let us write the normalized P* in dimensionless form: 



p* 



\ 



(35) 



/ 



ct/7 + /3/7 
a/P 

where Z = 1+ a/ft + ft/a + a/j + ft/ r y differs from Z in Equation (flUj) by a factor 
of aft'y. The stationary probability currents can be easily computed from here. There 
are only two independent ones, e.g., K*\ = a/Z and K*\ = P/Z. For these simple 
currents, there is only one associated graph, each with an irreversible loop, illustrating 



Equation (114p . Again, the details are in Appendix B 



Turning to the dynamic equivalence classes associated with this process, the 
simplest is to add an arbitrary A that satisfies Equation ( |27j) . with 6 free parameters 
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(denoted by e^...^ here): 

/-(ei + e a + C3)f e, e 2 ^ e 3 g \ 

A= £ if -(ei + e 4 + e 5 ) e 4 ^ e 5 £ 

£2 J £4 -(e 2 + e4 + e 6 )^ e 6 £ 

V e 3f e 5 e 6 ^ _( e3 + e5 + e6 )£ y 

Being the most general A, the new transition matrix W+ A does not provide transparent 
insight. Nevertheless, we can draw some interesting conclusions. For example, while the 
spectrum of W is complex for some a, (3, and 7, the spectrum of W + A will be real 
for sufficiently large e's. Details for the case with only e 4 7^ are found in Appendix 
B. While this modification corresponds to the rather innocuous addition of backwards 
hops for the particle, more drastic changes can be made, such as €3 which allows particle 
pair creation/annihilation transitions! We reemphasize that, in all cases, none of these 
modifications will affect the NESS distribution of probabilities or their currents. 

Finally, we illustrate the effect of A on entropy production. Since the original W 
does not allow "backwards" transitions, the production rate is infinite. The addition of 
A changes this result to 

= Z' 1 {(3 In (1 + /3/e a ) (1 + /3/e 2 ) + a In (1 + a/e 5 ) (1 + a/e 6 ) 

+ (a + P) In (1 + (a + (3) /e 4 )} . (36) 

Interestingly, this quantity is independent of €3, a result that follows from the lack of 
transitions between configurations 1 and 4 in the original dynamics. 



5.2. Models of mass transport 

Motivated by a wide range of physical problems, simple models of mass transport have 
been introduced. As extensions of TASEP, these models are defined by more general 
rates for moving masses from site to site. Recently, the exact steady-state distribution 
for a large class of such models was found, so that they serve as good illustrations of 
the framework presented above. Here, we will consider the most elementary model, the 
zero range process (ZRP) [28J: discrete (but otherwise arbitrary) masses transported 
around a ring of discrete sites. Specifically, let each site s (s = 1, L) of a periodic one- 
dimensional lattice be occupied by an integer valued mass m s so that a configuration, d 
or just i, is specified by the set {m s }. Thus, w\P* in Equation (fl2l) would assume 
the form P* ({m s }) w ({m s } — > {m^}). In the random sequential updating version, 
a transport event takes place during an infinitesimal time interval dt, at some site 
k: A portion, /x, is chipped off from and added to m k+ i, according to a given 
conditional rate q(fi\mk). In other words, the transitions w ({m s } — > {m' s }) connect 
only configurations with all m's being identical except for a single pair. To be explicit, 
we have, 

m s = m' s for s 7^ k, k + 1 (37) 

while 

m k -m' k = ji; m k+1 - m' k+1 = -/i (38) 
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occurs with rate q (n\m k ). This particular "element" in P* ({m s }) w ({m s } — > {m' s }) is 
just 

g(m fc - m' k \m k ) P* ({m s }) 5 (m k + m k+1 ,m' k + m' k+1 ) J~J 5(m s ,m' s ), 

s=ik,k+l 

where S (•, •) is the Kronecker delta. Clearly, the total mass M = m s is conserved 
and, with finite L as well, the configuration space is both discrete and finite. 

In general, finding P* from a given q can be quite difficult. Fortunately, it can be 
found for a wide class of such models [29]: If (and only if) q (fi\m) can be cast in the 
form 

q (ji\m) =g(p)f(m-fi)/f (m) , (39) 

where g, f are two arbitrary positive functions, then P* is factorizable. Explicitly, 
P*(W) = Z-'Usfims), where Z = £ {OTs} $ (M, £ s m s ) ]J S f (m s ) is a 
normalization factor. In configuration space (hypercubic L-dimensional lattice {m s }), 
the steady-state current 

K* ( W -> {m' s }) = P* ({m s }) w ({m s } -> {m' s }) - P* ({m' s }) w ({m'J -> {m s }) (40) 

can exist only on planes spanned by m k and m^+i, i.e., those specified by Equations 
( |37j) and ( 1381) . The expressions reduce considerably, since the /'s from g and P* cancel. 
If > m' k , then the current is easily understood as the result of moving \x > 0: 

K* (m k , m k+1 -> m' fe , m' fe+1 ) = Z" 1 ^ (/i) / (m fc - /i) / (m s ) 

s^k 

= Z~ l g (m k - m' k ) f (m' k ) JJ / (m s ) 

s^k 

Alternatively, we can keep P* implicit and write 

K* (m k , m k+1 -> m' k , m' k+1 ) = g (m k - m' k ) k P* ({m s }) . 

In either case, a simple mnemonic for K* is: "Replace / (m k ) by (7 (m k — m'^) / (m^.) in 
P*." For completeness, we should point out that, for m k < m' k , there is no transition 
across the (k, k + 1) bond, of course. However, the net current is still non-trivial, due 
to the difference of terms in Equation (I40p . Thus, we write 

K* (m k , m k+1 -> m k , m k+1 ) = - Z' l g (m k - m k ) f (m k ) JJ / (m'J 

s^k 

= -g{ m ' k -m k )t^P*{{m> s }). 
A more general version of this system, highlighting other subtleties of mass transport 



models, can be found in Appendix C 



Before ending this subsection, let us illustrate how Equation ([27j) can be exploited 
here, to define other transition rates which will lead to the same {P*,K*}. Since P* 
is explicitly known, we can easily add further mass transports events without altering 
the final NESS. For example, in addition to chipping fi from m k and moving it to m k+ \ 
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as above, we may chip off amount of mass, jl, from some other site, £, and move it 
to another site, £', with an arbitrary rate w, provided that we also make the reverse 
move with rate wf (me — jl) / (me' + jl) / f (me) f (m^), according to Equation ( 1271) . Of 
course, jl cannot be completely arbitrary, since it must lie within [0,m£]. Further, many 
masses can be chipped off at various sites and moved simultaneously (i.e., allowing moves 
away from planes in the L-dimensional space), so that almost arbitrary configuration 
changes can take place. The NESS, as we characterize it, will remain invariant provided 
the reverse moves occur with a rate obeying Equation ( 1271) . It is clear that such "long- 
ranged" moves will be much more efficient in Monte Carlo simulations of the true P*. 
This is especially critical for systems that display a condensate, since, on the one hand, 
P* is translationally invariant (or totally symmetric under permutations of s), while, 
on the other hand, typical configurations are associated with a condensate located at 
a particular site. With the original rates, the tunneling times for the condensate to 
"move" from one site to another, required to restore the invariance of the true P*, are 
prohibitively long J30j . In contrast, with arbitrarily long range jumps in configuration 
space, it is possible to design transition rates that readily sample the symmetric P*. In 
this sense, such transitions are comparable to the cluster algorithms which proved so 
useful in simulation studies of equilibrium systems |25j . 

5.3. Ising models coupled to two thermal baths 

Since its inception 80 years ago, the static Ising model has attracted the spotlight 
of equilibrium statistical physics. To simulate dynamic behavior when the system is 
coupled to a thermal reservoir, many researchers proposed an assortment of transition 
rates that, not surprisingly, obey detailed balance 'with respect to' the famous 
equilibrium distribution. The two most important classes are physically motivated: 
Glauber dynamics [31 J involving flips of individual spins and Kawasaki spin-exchange 
[32] which conserves the total magnetization (or particle number, in the lattice gas 
formulation). Given that this model is so well studied in equilibrium, it offers itself 
naturally to investigations of non-equilibrium steady states. Since the 80's, a large 
variety of rates that violate detailed balance have been introduced, leading to an 
extremely rich array of NESS-Ising models. The uniformly driven lattice gas [33j [T], 
based on Kawasaki dynamics, is a prime example. Another class of such NESS is 
the two-temperature Ising model. Based on Glauber spin-flip dynamics, it is coupled 
to thermal baths at two different temperatures [T9J [20j [T]. Needless to say, multiple 
baths can also be considered, but we focus on just two, for simplicity. Even with 
this restriction, the possible implementations of such couplings are seemingly endless 
(2TJ [TJ . For our purposes here, we consider a one-dimensional chain with spins coupled 
alternately to the two baths. The advantage of this model is that, in the steady state, all 
correlations are known [23] and energy fluxes through the system have been computed 
[2T] . Indeed, the full time-dependence is also accessible [31] ! 

As before, we denote the sites of a periodic one-dimensional lattice by s (s = 
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1,...,L). Each site carries an Ising spin a s , so that a configuration, C, or just i, is 
specified by the set {a s } (or just {a}). Of course, the energy of a configuration is given 
by the usual form: — J^2 S a s a s+ i. Choosing, for convenience, an even L, we couple the 
spins at the even/odd sites to two different thermal reservoirs, with temperatures 

T e = T even s ^ Todd s = T . (41) 

Using specifically heat bath dynamics, the master equation for this system can be written 

as 

dtP (M , t) = J> ({</} - {a}) P ({a'} ,t)-w ({a} - {a'}) P ({a} , t)] 
{*>} 



(42) 



where 



u>(K}-M) 




1 + IsOS 



a 



s-l 



s+l 



and 



7 S = tanh 



2J 

k B T s 



and wq just sets a scale. Note that this dynamics is quite special, since the transition 
rates w are independent of the state of the spin to be flipped. Furthermore, with random 
sequential updates, only a single spin is affected at a time. Thus, the probability currents 
exist only along the "edges" of an 2 L -dimensional cube (the corners of which represent 
the full configuration space). To proceed, it is useful to define a spin-flip operator JF S , 
acting on a spin configuration {a} : 

Fs W = {•••Cs-i, -a s ,a s+1 , ...} 
which just flips the spin at site s. Now, we can write the net stationary currents, e.g., 



K*{{a}^T s {a}) =^ 

_ mi 

2 



1 + ls°s 



0~ s -l + 0- s+ i 



s+l 



p* (W) 

p* {a}) . 



(43) 



In the limit T even s — >■ T Q dds, all the 7's are the same and P* ({cr}) oc 
exp [(J/k B T) J2 s a s a s+i], so that fT* = 0. For the non-equilibrium simple and 

compact form for P* ({cr}) is yet to be discovered. In particular, from the known 
correlations, we can deduce that the "effective Hamiltonian" would involve long range 
interactions [22]. Probably the most convenient representation, due to Hilhorst [35], is 



where 



1 + X L 



E 



n 



1 + Xt s t. 



s's+l 



II 1 - 



h, 



n 

.odd s 



1 + h r s a s 



(44) 



A = tanh 



tanh 1 V7e7o 



and h e>0 = J (j e + j a ) /2^ 0fi 



Probability currents as principal characteristics 



22 



Note that this representation requires performing a configuration- like sum over an 
auxiliary set of spins {r}, so that it resembles the partition function for a one- 
dimensional Ising chain in an inhomogeneous magnetic field. Thus, it is quite involved 
to find the probability currents K*. Not surprisingly, they are not zero in general. 
Deferring details to elsewhere, let us provide some illustrations here. 

If the neighboring spins are opposite (<7 s _i 7^ cr s +i), the transition rates, w, for 
flipping cr s are the same in both directions. For an equilibrium system, such a local 
environment implies P* ({er}) = P* {J-'s {o~}), so that the steady-state current is trivially 
zero. By contrast, we can show |36j that 

P*(W)^P* {F s {a}) 

for, e.g., the configuration {a} : o~ s -i = —1 with all other spins positive, when 
T even s 7^ T odd s . Thus, the net current K* ({a} — ► T s {a}) is non-zero. Of course, it is 
proportional to the difference of the two temperatures and vanishes in the equilibrium 
limit. Its explicit form is rather complex without being especially illuminating, and 
will not be displayed here [36J. Instead, let us present the sum of these currents, over 
all spins not involved, i.e., {a} = {..., cr s _ 2 , cr s+2 , ...}. For this, we define a reduced 
probability 

p(a- a _i,<7 aj <7 a+1 ) = 53P({<7}) ( 45 ) 
so that the sum 



is simply 



Wo 

2 



1 - IsPs 



JC(a x 



08-1 + 0. 



-cr a |c7 a _i,cr a+ i) = K* ({a} -> T s {a}) 
{*} 



s+l 



P{(T s -l,(T 8 ,(Ts+l) ~ Y 



1 + 7^. 



0~ s -l + CT< 



+ 1 



P {(7 s -i, -<J S , 



From (jHJ), it is easy to show that, in the L — > 00 limit (where A L — > 0), 



p (cr s _i, a s , a s+1/ 



1 + A/i s _ir s o- s _! 



1 + Xh s+1 r s a s+1 



[l + Xh e h (cr s _! + cr s+ i) cr s + (A/i s _i) 2 cr s _icx s+ i] 



where we have used h s -i = h s+ i and h s -\h s = h e h Q . As a result, 

w 



K 



[2A/i e /i - 7s (1 + A 2 /^)] K (a s _i + o- s+ i) /2} 



Note that the last factor assumes the values ±1,0 only. After some algebra, we find a 
particularly simple result: 



\K (a* 



— O's\0's-li &S+1, 



\le ~ 7o| ( wo/16) for cr s ._! = a s+ i 







for <r a _i 7^ 0- 



(46) 



s+l 



Note that this current sum vanishes for a s -\ 7^ c s +i, even though individual currents 
may be non-zero. Such a result is not surprising, given that the sum restores an 
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underlying up-down symmetry. Finally, we remark that all quantities for finite L are 
available, and this expression turns out to be exact for any L. Details of these findings, 
as well as extensions such as computing current loops, will be published elsewhere [36J. 

We end this subsection by illustrating how the currents can be used for computing 
fluxes. In particular, we recover the results of [21] concerning the energy flow into/out 
of the spin chain from the bath with the higher/lower temperature. When one spin (say, 
cr s ) is flipped, the energy change in the system is just 

AH (a s — > -a s ) = 2Ja s (a s -i + a s+1 ) . 

Since this change is independent of all the other spins, we can simply multiply it with 
/C (a s — > —a s ) to obtain the rate of net change (in the steady state). The result is 
wq J (7s-i — 7s) (o~ s -i + °~s+i) 2 /8, so that the configurational average is wqJ (7^-1 — 7 S ). 
Note that, if T s > T s _ 1; then 7 s _x > 7 S so that there is an average energy flow "into" 
the system due to flipping the spin at site s. Thus, the net flux through our system is 
w J |7 s _i — 7g I /2, a result phrased as "the energy flow from an even to an odd site" 
in [21]. Interestingly, these energy fluxes are directly related to the entropy production 

[221 EE]. 



5.4- NESS Gaussian distributions 

Another exactly solvable system is M "particles" (or simply, degrees of freedom), 
subjected to linear forces (e.g., generalized coupled harmonic oscillators) and Gaussian 
noise. Although the configuration space of this system is continuous, it is sufficiently 
simple that both P* and K* are explicitly known, serving as an elegant illustration of 
our formalism. It can be regarded as the continuum limit of a particular biased random 
walk on, say, a hypercubic lattice. Specifically, the jump rates have a bias which depends 
linearly on the location of the walker from the origin and may be anisotropic. 

We begin with a Langevin equation for the degrees of freedom ( "coordinates" ) £ M 
(fj, = 1, ...,J\f) and the associated Fokker-Planck equation for P tj . The latter is the 
master equation for this system, with a stationary distribution known to be Gaussian. 
If detailed balance is violated, the stationary probability currents can nevertheless 
be computed. This subsection will be devoted to the highlights, with details left to 



Appendix D In this spirit, we will use a compact notation of vectors (e.g., £) and 



matrices (blackboard bold font) here, while in Appendix D all indices will be explicitly 
displayed. 

Consider a Langevin equation with linear deterministic forces: 

d t £(t) = -F£(t)+ff(t) (47) 

and Gaussian noise (uncorrelated in time): 

(tf(t)> = 

(ff{t)®rj{t')) =N5{t-t'). 
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Here, F is an arbitrary real matrix, except that, for stability and to model dissipation, 
the real part of its spectrum must be positive. We denote the eigenvalues and the right 
and left eigenvectors of F by 

Fm/ = Aj«j and tf/F = Aj?Tj (48) 

with Re A/ > 0, and I = 1,2, ...,JV. Meanwhile, the noise matrix N must clearly be 
real symmetric, as well as positive. For example, if the discrete random walker has 
anisotropic jump rates, then N is diagonal but not proportional to the unit matrix. 
The associated Fokker-Planck equation reads 

d t P (f, t) = V • 1 - V + wA P (£ t) (49) 

with the probability currents 

£ = _ |5v + F|jp(£t). (50) 

If the matrices N and F are constrained such that N _1 F is symmetric, then we can 
define 

V = 2jfc B TN _1 F (51) 

and see that this is just a system of coupled simple harmonic oscillators, subjected to 
(a general Gaussian) noise. Such a system will settle into thermal equilibrium, with the 
familiar Boltzmann distribution: 



P eq (t) oc exp 



2k B T 



as well as a trivially vanishing K. Of course, Equation (I5T1) with symmetric V is just the 
detailed balance condition here and expresses the usual fluctuation-dissipation relation. 

However, if we insist on studying the general case where N _1 F is not symmetric, 
we will encounter a NESS with non-trivial K's. (For M = 3, we would write such a 
K as curl of some vector field.) To see this explicitly, we recognize that the stationary 
distribution is still a Gaussian Wt 



P* (e) ex exp 



(52) 



2 

where C is the correlation matrix: 

'£®£) = C. (53) 



To find C in terms of N and F is straightforward, but not trivial (Appendix D). A 
convenient expression is: 



A/ + Aj 

which is manifestly symmetric. 
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Turning to the stationary currents, we see that the second term in 

= _ J5v +F |jp*Q (55) 

can be written as -FCVP*. Thus, K* = {FC - N/2} VP*. To see that FC - N/2 is 
indeed antisymmetric (so that V-K* = 0), we only need to recall that uj are eigenvectors 
of F, so that 

V ' + J 
The final result is quite simple: 

K* = A VP* , (56) 



where 



1 y- h_ 

9 2^ \ r 



AJ (v i Nvj)u i ®uj (57) 



2^Xi + Xj 

is manifestly antisymmetric. Of course, we can evaluate the gradient in Equation (156 
explicitly, with the result 

K* = - ( AC^tj P* . (58 



One interesting consequence is that the expression in Equation (ITTj) can be neatly 
evaluated: 

2 1 



^Tr {AC' 1 A T } 



Lastly, let us consider "generalized detailed balance" : What is the class of dynamics 
that will lead us to a given NESS, i.e., {P*, K*} or more specifically in this case, {C, A}? 
Since the configuration space here is not finite, it does not seem trivial to provide the 
entire class of allowable dynamics. However, since we know P* explicitly, we can exploit 
Equation (J27J) and add / d£' j A(£ £')P{£',t) - A(|\ <f)P(£, t)} to the right hand side 

of Equation (|49l) . where A is any function that is non- negative and satisfies 



A(e,e 



exp 

Unfortunately, this expression is too general to be illuminating. To provide some 
insight, let us illustrate this freedom by studying a small subset of the allowable 
modifications. Let us frame this question as the "inverse" of the usual one. Ordinarily, 
we are given the matrices {F, N}, and are asked to find the NESS (in this case, {C, A}). 
The inverse question is: Given {C, A}, what can be said about {F, N}, i.e., what 
dynamics will insure that we arrive at the above {P*, K*}7 Clearly, the key equation is 
just FC — N/2 = A. Specifically, we are allowed to choose any N, provided it is a valid 
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noise correlation (real, symmetric, positive matrix) and construct the dissipative forces 
F according to 

F = CT 1 {N/2 + A} . (59) 

It may also appears as if we can choose any matrix F and simply compute N from 

N = 2 {FC - A} . (60) 

A closer examination reveals a difficulty with this naive approach. For fixed {C,A}, 
not every choice of F will lead to a symmetric N (a necessary condition for the noise 
correlation) ! This issue is related to the discussion of Equation ( 15TI) . 

We end this subsection with two comments. The concerned reader may ask: Why 
can we choose an arbitrary N but not any F? The difference lies in that N has fewer 
parameters than F. There is an intimate connection to the general case where only the 
symmetric part of Wf can be chosen freely. Since N must be symmetric, it evidently 
contains all the freedom of choice available (in this subclass of dynamics). By contrast, F 
appears to have no particular symmetry; yet the antisymmetric part of CF is completely 
fixed (by A). Finally, to make contact with ordinary detailed balance, it is easier to 
consider Equation (1601) and follow how it reduces for the textbook case of an equilibrium 
system with the most trivial dynamics: We have (i) a Hamiltonian Ti = £V£/2; (ii) forces 
k ^— V7Y^ where k is a dissipative coefficient; and (iii) C = (V / IcbT)^ 1 . Thus, A = 0; 
and N = 2FC is just 2KksT x the unit matrix. 



6. Summary and Outlook 

In this paper, we propose that a non- equilibrium steady state be characterized not only 
by P* (C), the stationary distribution for finding the system in configuration C, but 
also K* (C,C), the net probability current from C to C. The pair {P*,K*} forms a 
natural generalization of having P eq as the quantity that characterizes all details of 
an equilibrium state. Within the context of a master equation approach, in which the 
dynamics of a system is defined by a set of transition rates w (C — > C), we recalled a 
graphic method for constructing the pair {P*, K*}, as well as the intimate relationship 
between K* and irreversible cycles associated with the rates. Assuming ergodic rates, 
this steady state is unique. 

Exploring the consequences of our postulate, we find that the converse - finding 
a set of w's which leads to a given {P*,K*} - is far from unique: The rates are only 
bound by a certain set of constraints. This condition is well known, to those who study 
equilibrium systems by computer simulations, as "detailed balance:" 

w (C -> C) P eq (C) - w (C -> C) P eq (CO = 

The proposed generalization to non-equilibrium systems is almost intuitive, namely, 

w (C -> C) P* (C) -w(C -> C) P* (CO = K* (C, CO • 

Since the only difference between this and the equilibrium condition is the value 
of the right hand side, we are "just as free" in choosing rates in either case. In 
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other words, the equivalence classes of rates are the same, except for one subtlety 
explained in Section HI A unified way of making this statement is the following. Two 
sets of rates belong to the same class provided their differences, A (C — > C), satisfy 
A (C -> C) P* (C) = A (C — > C) P* (C). A further consequence of these degrees of 
freedom is that the entropy production (associated with the reservoirs coupled to our 
system, driving it to a NESS [13]) can be made infinitesimally small or arbitrarily large. 
We also proposed that S (K*) 2 can be exploited to measure of how "far" a NESS is 
away from equilibrium. Of course, it would be better if a dimensionless quantity could 
be identified, so that common phrases like "systems far from equilibrium" can be given 
a universal and quantitative meaning. Finally, a number of examples were presented, 
illustrating these ideas. 

While we hope that this work provides a fresh perspective on non-equilibrium steady 
states, we are aware of much room for improvements and further investigations. Let us 
end with a brief outlook. 

For simplicity, we focused mainly on a master equation with continuous time and 
discrete configuration space. Generalizations to systems with continuous configuration 
space should be possible, though there may be non-trivial mathematical obstacles. We 
believe that the example in 15.41 represents a possible starting point. The more general 
case is the Fokker-Planck equation, which admits transitions between configurations 
{4>(x)} that are infinitesimally close: 5(f) (x). Much work has been devoted to this 
equation for NESS [38], and the resulting conclusions should be exploited. On a 
different note, formulations involving discrete time present a different challenge. Since 
the transition "rates" are now conditional probabilities, their proper normalization will 
lead to further constraints on the choice of w's, beyond those listed in Section HJ) . 
However, we anticipate this complication to be minor. 

As more serious question, yet to be answered for arbitrary NESS's, is the existence 
and uniqueness of thermodynamic limits. For example, there is a belief that, for the 
driven Ising lattice gas in two dimensions, different steady states will be reached in 
this limit, depending on the aspect ratio of the system [39]. Perhaps devoting some 
attention to the distribution of probability currents will facilitate this quest. This topic 
is also intimately related to the issue of the "micro-macro connection." Starting as a 
vague notion of "coarse graining," this connection has gained much substance through 
the renormalization group (RG), especially in the study of critical phenomena. For 
systems in equilibrium, the distribution (P eq ) is uniquely linked to a Hamiltonian (7i) 
by the Boltzmann factor, and the progression from a microscopic, via a mesoscopic, to 
a macroscopic description can be cast in the language of RG flows in the space of TVs 
or P's. We believe that, for NESS, the flow in the space of K*'s will play an equally 
important role, on a par with the flow of P*. To illustrate this thought, let us highlight 
two seemingly similar NESS systems in which the RG flows end at very different types 
of fixed points. Specifically, for the uniformly driven Ising lattice gas [33] , the RG fixed 
point is a "genuinely non-equilibrium" dynamic functional [10] which cannot be written 
in the form of any equilibrium system with an Ti. By contrast, if the same system is 
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driven randomly, the fixed point corresponds to an equilibrium uniaxial system with 
dipolar interactions [41j . This example suggests that the RG flow of K* leads to some 
non-trivial fixed point in one case and, in the other case, ends at the "trivial" fixed 
point K* = 0. One hint from physics is that there are global (particle) currents in the 
first system, which cannot vanish under the RG. In contrast, non-trivial microscopic 
probability currents in the second system are most likely local, so that they vanish 
under RG transformations. Work is in progress to provide solid foundations for these 
promising notions. 

Entropy production is another venerable issue. Through its coupling to more than 
one energy reservoir, our system can be regarded as an agent which redistributes energy 
from one or more of these reservoirs to the rest. The most intuitive and simple example 
is given in Section 15.31 above - an Ising model coupled to two thermals baths at different 
temperatures. Heat flows from the hotter bath into the spin system and then to the 
cooler bath. In the steady state, all thermodynamic quantities of the Ising system are 
constant on the average, including this heat through-flux. Thus, there is a constant rate 
of energy redistribution associated with "the environment" (or "the medium"), a term 
we use to refer to the totality of all the external reservoirs. Associated with this energy 
redistribution should be an entropy production of the environment. In general, given the 
detailed dynamics of the environment and how it couples to our system, it is possible, 
in principle, to compute this rate of entropy production. However, in our framework, 
we start with a master equation. Expressed in terms of transition rates, it captures the 
dynamics of the system and how the system couples to the environment, but carries 
no information about the dynamics of the latter. As a result, it is not clear how we 
can unambiguously find the rate of entropy production of the medium. Nevertheless, 
it seems possible to define an entropy production, associated with the system-medium 
coupling. Relying on the context of chemical reactions, Schnakenberg proposed such 
a definition, which depends only on the io's in the master equation. Thus, it is a 
natural candidate for us to exploit here (see Section H~2l) . even though its relation to the 
true entropy production of the medium is not definitive. There may be other, better 
grounded, definitions [12]. We believe it will be worthwhile exploring them also in the 
context of K*. In a related vein, we mention fluctuation-dissipation and work-energy 
theorems associated with NESS [43] . Apart from considerable theoretical interests over 
the last two decades, recent experimental measurements in a variety of systems have 
contributed to the excitement focused on these issues |Sj. Whether the probability 
currents provide novel insight here deserves to be studied. 

Finally, let us return to the analogy with electrodynamics, in which a parallel is 
drawn between electrostatics/magnetostatics and equilibrium/non-equilibrium steady 
states. Of course, we are aware of considerable conceptual difficulties, especially for the 
general case where transitions between all configurations are allowed. However, for cases 
with a local dynamics (e.g., Glauber spin flips for an Ising model or a Fokker-Planck 
equation in continuous configuration space), an induced metric can be defined jl5] so 
that the probability current is just a vector field. Since K* is divergence- free, it can 
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be expressed as a (generalized) curl of a vector potential and 'magnetic fields' can also 
be defined [15]. It behooves us to ask whether such fields might be meaningful, and if 
their properties can be exploited. Pursuing this parallel with electrodynamics far into 
the future, we may seek an underlying gauge theory [16] and see if it might provide a 
unified description of all time-dependent phenomena in statistical mechanics. 

The comments above highlight just a few of the issues surrounding the central quest: 
providing a sound framework for the statistical mechanics of non-equilibrium steady 
states, on a par with the Boltzmann-Gibbs approach for equilibrium systems. We hope 
that our proposal - to bring the microscopic probability current distribution, K*, onto 
center stage, as an equal partner with the microscopic distribution of probabilities, P* 
- will generate further explorations and discussions in the pursuit of our quest. 
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Appendix A. Kirchhoff's problem and its relation to NESS 

In most expositions of the graphic approach to finding the stationary distribution from 
the master equation, Kirchhoff's solution [15] to the electrical network problem is 
referenced. Here, we offer some comparisons and contrasts. 

At first glance, it appears that both problems have the same underlying network 
structure (N vertices and N (N — 1) /2 edges) and the same number of parameters: 
2 x N (N — 1) /2. For electrical networks, Kirchhoff [15] posed a resistance (Rij) and 
an electromotive force (£ij) between the nodes % and j. The former/latter are strictly 
symmetric/antisymmetric. Naively, since we have w\ (with i ^ j), we may think of 
identifying the symmetric and antisymmetric parts of w with R and £ , and perhaps 
formulate an exact mapping between the two problems. 

On closer examination, differences emerge. To seek the NESS, we need to solve 
N equations for N unknowns (P*). No other equations need be solved to obtain the 
N (N — 1) /2 currents; they are just linear combinations of the P*'s. By contrast, by 
focusing on the currents (iy, also antisymmetric) as the unknowns, Kirchhoff needed 
far more equations. In his paper, there are N equations for the nodes: 







3 

for each i, and N (N 

£k\k2 



3) /2 + 1 equations for the loops: 

+ £k n k x — Rkik 2 hik2--- ~ Rknki^knki = • 
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Of course, we recognize that the node equations correspond to the master equation for 
the steady state and that loop equations allow us to define a single-valued function 
associated with each node. Traditionally, this is called the "potential," Vi, which can 
be defined through the differences 

V{ Vj = £{j Rijlij . 

Clearly this is not the only possible definition of the single-valued function, since any 
well-behaved function of Vi will produce another single-valued function. Even if we 
demand linearity between V and J, the above is unique only up to two parameters: a 
constant shift and an overall scale: V — > A (V — v). If we were to try to identify P* 
with the potential, we will need the shift to insure positivity (since potentials can easily 
be negative) and the scale for normalization of P*. But the shift is arbitrary (apart 
from a lower bound), so that the ratios (Vi—v)/ (Vj — v) are not fixed. In this sense, 
there is no way to generate a unique quantity like P* from the currents of the Kirchhoff 
problem. 

Forcing the issue from the other direction seems more feasible, since there is no 
problem choosing the potentials Vi to be P*. From the definition K*\ = w\Pj — WjP*, 
we can write 

w- I 

p* _ _1 p* j 

i w) j w) 1 

Formally, we can substitute the expression in Equation ( flOj) for P* and write the first 
term ( w\P* /m* ) as G\ \{w)\- Then, we have 



P*-P*= [Gj-G}] 



1 1 



Taken at face value, we can exploit this equation to identify 
P* -V, 

Li i — Lr j — > tij 
"1 _> 

The lack of a one-to-one mapping between {R, 8} and {w} can now be summarized 
succinctly. From the former, we can find a unique solution for the currents, but not the 
charges (or potentials). Even if we could and these were identified with K* and P*, we 
know that there is a whole class of w's which leads to the very same {P*, K*}. Finally, 
in passing, we mention that there are intimate connections between Kirchhoff 's problem 
and the percolation problem. Details may be found in the context of an extensive review 
on the Potts model [37]. 
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Appendix B. A maximally asymmetric iV 
sites 



4 example: TASEP with two 



For a totally asymmetric exclusion process (TASEP ) on a one-dimensional chain with 
open boundaries [TT], [18] with two sites, we show all the configurations in Figure HI 
From the transition matrix, Equation ( 1341) . we see that there are only 5 non-zero arrows 
in the full space, shown in Figure [5](a). As a result, there are very few non-trivial 
directed trees. In Figure EJ^b), we show the only tree directed towards vertex 1 that is 
non-zero. For vertex 3, there are two trees, shown in Figure [5]^c) and (d). Thus, we 
arrive at, respectively, P{ 2 3 4 oc PjP, atflj, aa(3 + /3/3a, aja , and Z in Equation ( TTOT) is 
just a 2 (3 + P 2 a + (a 2 + a (3 + /3 2 ) 7. To illustrate Equation (THj) . we consider, e.g., K*\ 
which is j3/Z. In graphic terms, we add a down arrow from vertex 1 to 3 to Figure E^b), 
forming an irreversible loop associated with II [C] = a[3j. Meanwhile, there is only one 
side branch, so that R — (3 here. Thus, K*\ = UR/Z = a/37/3 / '(aPjZ), in agreement 
with the above. 

From these X*'s, we can evaluate the particle current flowing into the system. 
Referring to Equation ( |20l) . it is clear that we just need to consider (1 — n\) n\, which is 
zero for all configuration pairs, except being unity for the (1,3) and (2,4) pairs. Thus, 

(J*) = K*l + K*l = (a + (3) /Z . 

In steady state, it is clear that this expression will also equal the current flowing between 
the two sites: K*%, or out of the system: K*\ + K*j. 

Next, we turn to the dynamic equivalence classes associated with this W. We first 
compute the matrix W, and then extract its symmetric and its antisymmetric part, 
resulting in: 
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(B.l) 



(B.2) 



We clearly recognize the values of the K*'s in the elements of A. 

Since S is a 4 x 4 symmetric matrix with constrained diagonal elements, we have 6 
degrees of freedom. Selecting just one of these, for the purposes of illustration, we can 
define a new matrix, S, by, e.g., modifying the (3,2) and (2,3) elements by an amount 
2e, and preserving the column sums 

/ -2p P P \ 

1 P -2(a + P)-2e a + p + 2e a 

p a + p + 2e -2{a + P)-2e a 

a a —2a 
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the new rate matrix W is given by 
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Comparing to the original matrix W, we can see that the new feature is a transition 
from configuration 2 to 3 which was absent in the original model. In other words, the 
particle is now allowed to jump backwards. The matrix A is nonzero only in the 2x2 
submatrix involving configurations 2 and 3: 



A = W — W — e 
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It is easy to check that it satisfies Equation (127)) . 

Let us now recalculate the particle current from the first to the second site, keeping 
in mind that the particle can now also jump backwards. However, the probability current 
between configurations 2 and 3 is unchanged, as one can see from a quick explicit check: 

7 \{a + P) 
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a + P J 7 

Since each term in the equation above is associated with a single particle hop, it is 
immediately obvious that (J*) also remains invariant, under the transformation from 
W to W. 

Allowing for backwards hops of the particle may appear rather innocuous. A more 
drastic modification would be to allow transitions between configurations 1 and 4 which 
are completely absent in the original model. This leads to 
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Again, the difference between the two sets of rates satisfies detailed balance with respect 
to P*. As a result, the probability current K*\, associated with transitions between 
configurations 1 and 4, remains zero. This also ensures that the physical current, (J*), 
remains unaffected. 

Since all of these rates involve at least one uni-directional transition, they are all 
associated with infinite entropy production. In order to obtain a finite value for the 
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entropy production, we need to start from a rate matrix where each forward transition 
is associated with a backward one. Starting from the original W and its associated S, 
Equation ( IB.lj) . and recalling that we have 6 degrees of freedom, parameterized by e%, 
e 2 , e 6 , the most general S would be 
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The associated W is given by 
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Evaluating the entropy production for this general form, we arrive at 
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which is clearly finite now. 



Appendix C. Currents in generalized mass transport models 

Here, we provide a further example associated with mass transport models, namely, 
the generalization from discrete to continuous masses, with parallel updates in discrete 
time. Thus, we consider a ring of discrete sites occupied by continuous m s 's. The 
evolution proceeds in discrete time steps, with /x s being chipped off from m s and added 
to m s+ i, for all sites at once. To distinguish the conditional probability here from the 
case in the main text, we denote it by here. Though the fundamental issues for 

this non-equilibrium process and the one with continuous time and random sequential 
updates are the same, there are a few differences, on which we will comment here. 

As in the previous case, P* is known for a wide class of such models[29j : If <ft {tA m ) ls 
of the form v (//) u(m — //)/ / (m) (where v and u are arbitrary non-negative functions, 
and / (m) = f v (/i) u{m — fi) dfi), then P* ({m s }) = Z^ 1 Y[ s f (rn s ), as in Section 15.21 
Again, wjP* in Equation (JT2l) will assume the form P* ({m s }) w ({^ s } —>■ { m ' s })- Unlike 
the random sequential case, 

w ({m s } -> {m' s }) = JJ (jJi a \m e ) 



Probability currents as principal characteristics 



34 



so that all the f's cancel. Thus, we are left with a simple product: 
P* ({m s }) w ({m s } - {m' s }) = Z' 1 J] v (//.) it (<r fl ) 

s 

where fi s (the mass to be moved) and cr s (the remaining mass) are fixed by 
m s = fj, s + a s 

m' s = /ig-i + a s (C.l) 

Thus, 

ZK* ({m s } -> {m'J) = JJ v (ji.) u (a s ) -J[v (ft s ) u (a s ) , (C.2) 

s s 

where ft s and a s satisfy m' s = ft s + a s and m s = p, 8 _i + a s . Let us end with a few 
remarks. 

First, it is clear that many different functions v and u will lead to the same 

convolution / = v * u, so that all of them are associated with the same P*. However, 

it is also clear that typically, K* will not be the same. Thus, our framework provides a 

natural way to distinguish these different NESS's, all of which share the same stationary 
p* 

Next, we should emphasize that this rather compact notation belies a subtlety: 
Starting with a pair of configurations {m s ,m' s }, it is not possible in general to find a 
unique set {/i s ,cr s }. First, Equation (IC.lj) . rewritten here as 

fJLs-l +c s = m' s 
fi s + cr s = m s 

is just a simple linear equation of the form 

/ 1 ... 1 \ / <ri \ / m[ \ 
1 1 ... \i\ m\ 



... 1 a L m' L 
\0 ... 1 l/\/iJ \m L J 

Now, the matrix on the left has a zero eigenvalue and cannot, strictly, be inverted. 
However, the zero eigenvector (+1/ — 1 for even/odd elements) is associated with mass 
conservation. A physical set of {m s , m' s } obeys Yl s ( m s ~ m ' s ) = an d has no component 
along this eigenvector. Thus, the existence of a solution is not in doubt. 

In general, the solution is not unique, however. With parallel update, the 
transported masses {/i s ,cr s } are determined by configurational masses {m s ,m' s } only 
up to an overall parameter /2. This extra "freedom" is associated with a simple physical 
process: The same configurational change can be achieved if we further chip off an 
identical amount, /2, from all sites and moved to the next, i.e., replacing {/x s ,cr s } by 
{fi s + p,, a s — p,}. The range of allowed p, is clear: None of the resultant /i's and a's 
can be negative. Thus, if we start with any "minimal" set of /Vs (i-e., at least one 
fi is zero), then the other sets can be generated by adding a positive ft, provided it is 
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no larger than the smallest m s in the system. In other words, from any valid solution 
(i.e., non-negative {fi s ,cr s }) of Equation fIC.lj) . /i can decreased (or increased) until the 
smallest fi s (or o~ s ) vanishes. In this sense, we should integrate over this allowed range 
of Ji when computing the first term in the current in Equation (1C.2j) . 

Once a {/x s ,cr s } solution is found, there is a simple way to construct the 
"complementary" set {fi s ,a s }, which satisfies 

jx s -i + a s = m s 
fx s + a s = m' s 

Indeed, we can access the minimal set immediately. Let fx be the largest of the /i's: 
jj, = sup {/i s } . 

Then, 

L~l s = (l- Ll s 

OS = OS + fU s + - LI. 

Note that the one parameter family of allowed {fl s ,a s } is not necessarily the same as the 
one for {/z s ,cr s }. As a result, a different integral is needed here to compute the second 
term of K* in Equation (1C.2I) . 



Appendix D. A simple example with continuous configuration space. 

In the following, we provide the detailed analysis of the Langevin equation with linear 
deterministic forces (generalization of a system of coupled harmonic oscillators). The 
configuration space is Af- dimensional space: with 11 = l,...,Af. The Langevin 

equation for (t) reads: 

where F{? is an arbitrary real matrix, except that, for stability, the real part its spectrum 
must be positive (dissipative forces). The noise r] is Gaussian with zero mean, i.e., 

fa"(*)> =o 
(r]^(t)r] v (t')) = N^5(t-t') 

where N^ v is real symmetric and positive. 

Let the eigenvalues and right/left eigenvectors of F be given by 

J>>? = A(/)<; 5>JF£ = A (I) i£ 

V (J, 

I =l,...,Af 

with 

Re A (I) > 0. 

We choose eigenvectors so they form a complete bi-orthonormal set, i.e., 

i ii 
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The Fokker-Planck equation is 

d t P (£,t) = {\d^d v + d^A P [it) , (D.l) 

fl,V 

where 

a = A 

differentiates all factors to its right. 

The stationary distribution is also a Gaussian 



P* (tj = W N ' 2 (detG)- 1/2 exp -i J^G^ 



(D.2) 



where 

g = cr 1 

with C being the correlation matrix: 

The explicit form of C is found in the standard way, except that matrices are involved. 
Thus, it is formally 

. du ( 1 \ / 1 
C= / — N 



2tt \iuj + ¥ J \iu + F y 

In the frame where F is diagonal, this integral is trivial. The result is, in terms of the 
explicit matrix elements: 

r „_sr u i nIJu j 

A,A(I) + A(J) ' 

where 

n ij = v i N>MUv u 

is just N in the new frame. 

The probability currents are displayed in Equation (ID.lj) : 

v ) 

so that the stationary currents are 

k** (i) = -J2 {\ Nia,d » + F ^ v \ p* (i 

To see that this is the (generalized) curl of a vector field, we exploit the Gaussian 
property of P* and substitute 

£ 7p * = _J2 C^ u d u P* 



Probability currents as principal characteristics ... 37 
into the second term. Thus, 

k*» (i) = \ F i° lu - \ N * V \ d » p * (£) ■ 

v \ 7 J 

Next, we want to show that the matrix in {...} is indeed antisymmetric. Recall that u 
is an eigenvector of F, so that 

For N^, we cast it in terms of N IJ , i.e., 

so that the terms in {...} combine easily. The result is 

V 

with 

A(i) - A(J) 7J 

"2^A(/) + A(J)^ iV 

which is manifestly antisymmetric. Note that this should not be confused with the 
A\ in Section HI 
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